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Abstract  -  Raman  spectroscopy  is  a  powerful  and 
effective  technique  for  analyzing  and  identifying  the 
chemical  composition  of  a  substance.  Two  types  of  Raman 
spectra  estimation  algorithms  exist:  supervised  and 
unsupervised.  In  this  paper,  we  perform  a  comparative 
analysis  of  five  supervised  algorithms  for  estimating 
Raman  spectra.  We  describe  a  realistic  measurement 
model  for  a  dispersive  Raman  measurement  device  and 
observe  that  the  measurement  error  variances  vary’ 
significantly  with  bin  index.  Monte  Carlo  analyses  with 
simulated  measurements  are  used  to  calculate  the  bias, 
root  mean  square  error,  and  computational  time  for  each 
algorithm.  Our  analyses  show  that  it  is  important  to  use 
correct  measurement  weights  and  enforce  the  nonnegative 
constraint  in  parameter  estimation. 
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1  Introduction 

The  Raman  effect  or  Raman  scattering  represents  the 
inelastic  quantum  scattering  of  a  photon  by  molecules  in 
liquids,  gases,  or  solids  [2-3].  When  light  is  incident  on  a 
molecule,  most  photons  are  scattered  elastically  so  that  the 
energy  or  frequency  of  the  scattered  photon  is  the  same  as 
that  of  the  incident  photon.  This  is  known  as  the  Rayleigh 
scattering.  A  small  fraction  (about  one  in  a  million)  is 
scattered  inelastically,  causing  the  frequency  of  the 
scattered  photon  to  be  different  (usually  lower)  from  the 
frequency  of  the  incident  photon.  This  is  known  as  Raman 
scattering.  The  frequency  change  is  due  to  the  change  in 
energy  levels  of  the  vibrational  or  rotational  energy  of  the 
molecule.  Therefore,  Raman  spectroscopy  is  a  powerful 
tool  for  analyzing  the  chemical  composition  of  liquids, 
gases,  or  solids  using  a  laser  [2],  [13-15],  [12]. 


A  Raman  spectrum  is  a  plot  of  the  intensity  of  the 
scattered  photon  as  a  function  of  frequency  shift.  The 
measured  Raman  spectrum  can  be  used  as  a  fingerprint  to 
uniquely  identify  the  chemical  composition  of  a 
substance.  Application  of  Raman  spectroscopy  to  analyze 
chemical  compositions  of  various  substances  has  seen  a 
rapid  growth  in  recent  years  [2],  [13-15],  [12].  This  is 
primarily  due  to  the  development  of  inexpensive  and 
effective  lasers  and  charge-coupled  device  (CCD) 
detectors  [2].  Raman  spectroscopy  is  also  popular 
because  measurement  collection  is  fast  and  does  not 
require  contact  with  the  chemical  substance. 

Suppose  we  have  the  measured  Raman  spectrum  of  a 
substance  and  we  are  interested  in  determining  the 
chemical  composition  of  the  substance.  The  measured 
spectrum  contains  various  error  sources.  Therefore,  it  is 
necessary  to  use  a  statistical  measurement  model  that 
expresses  the  measurement  as  a  function  of  the  true 
spectrum  and  dominant  error  sources. 

Raman  spectrum  estimation  algorithms  are  of  two  types, 
supervised  and  unsupervised  machine  learning  algorithms 
[11].  In  the  supervised  approach,  a  library  of  reference 
Raman  spectra  are  used  and  the  true  target  spectrum  is 
expressed  as  a  linear  combination  of  the  reference  spectra. 
Each  reference  spectrum  is  assumed  to  be  error-free.  In 
practice,  this  is  not  feasible.  If  the  errors  in  a  measured 
reference  spectrum  are  very  small  compared  with  the 
signal  values,  then  it  is  a  good  approximation  to  treat  the 
measured  reference  spectrum  as  error-free.  Otherwise,  one 
must  model  the  errors  in  the  reference  spectra.  Supervised 
algorithms  assume  that  the  library  contains  all  reference 
spectra  that  may  be  encountered  in  data  collection.  A 
supervised  algorithm  estimates  the  nonnegative  expansion 
coefficients  or  mixing  coefficients  using  the  reference 
spectra  and  a  statistical  measurement  model.  The 
unsupervised  approach  estimates  the  spectra  and  mixing 
coefficients  directly  from  measurements. 
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This  paper  examines  estimation  of  Raman  spectra  using 
the  supervised  approach  and  performs  a  comparative 
analysis  of  five  Raman  spectra  estimation  algorithms: 

(i)  classical  weighted  least  squares  (CWLS), 

(ii)  nonnegative  weighted  least  squares  (NNWLS), 

(iii)  fast  combinatorial  NNWLS  (FCNNWLS), 

(iv)  block  pivoting  NNWLS  (BPNNWLS),  and 

(v)  NNLS  or  NNWLS  using  generalized  likelihood  ratio 
test  (GLRT). 

We  use  simulated  data  and  perform  Monte  Carlo 
simulations  to  calculate  measures  of  performance  (MoP) 
for  each  algorithm.  MoP  used  in  this  study  include  bias  in 
the  estimator,  root  mean  square  error  (RMSE),  bias  in  the 
measurement  residual,  RMSE  for  the  measurement 
residual,  and  computation  time. 

These  algorithms  are  implemented  in  a  software 
simulation  benchmark  system  designed  specifically  for 
analysis  of  detection  algorithms.  Each  algorithm  is 
inserted  into  the  benchmark  software  using  a  well  defined 
interface.  The  benchmark’s  native  language  is 
MATLAB 1 .  For  a  fair  comparison  of  run  times,  the 
algorithms  were  all  coded  in  MATLAB".  Nothing 
precludes  an  algorithm  from  being  implemented  using  a 
language  that  can  be  incorporated  into  the  MATLAB" 
environment,  i.e.,  Java,  C/C++,  or  FORTRAN.  However, 
for  purposes  of  comparison,  MATLAB "  was  used  for  all 
algorithms. 

The  outline  of  the  paper  is  as  follows.  Sections  2  and  3 
describe  the  measurement  model  and  measurement 
function  for  Raman  spectra,  respectively.  We  summarize 
various  Raman  spectra  estimation  algorithms  in  Section  4. 
Finally,  Sections  5  and  6  present  numerical  results  and 
conclusions. 

2  Measurement  Model  for  Raman 
Spectrum 


s  b 

where  ni  ,  nj  ,  and  g  represent  the  signal,  background 
noise,  and  Gaussian  noise,  respectively.  The  variables  tl] 
and  il-  are  modeled  as  discrete  random  variables  (RVs) 

s  b 

whereas  g  is  a  continuous  RV.  We  assume  that  lli  ,  H-  , 
and  g.  are  independent.  The  noise  g.  is  introduced  by  the 
on-chip  amplifier  and  is  modeled  as  Gaussian  with  mean 
m  and  variance  cr2 

gi  ~  N(gi;m,<J2),  (3) 

E{(g,  -  m){gj  -  m )}  =  Syff2.  (4) 

The  background  noise  It’  is  Poisson  distributed 

«,b  ~^Poisson  («/b;4b),  (5) 

where  /f’  represents  the  expected  number  of  counts  for 

the  background  noise 

E{nb)  =  4\  (6) 

and 

Ppoisson(X^)  =  ^^’  x  =  0,1,2,...  (7) 

x! 

We  note  that  the  variance  of  the  Poisson  distribution 
Ppoisson  (x’  's  also  A-  The  signal  is  Poisson 

distributed  with  parameter  A} 

Poisson  (8) 

where  A’  represents  the  expected  number  of  counts  for 

s  b 

the  signal.  Since  n]  and  ll  are  assumed  to  be 

s  b 

independent,  ni  +  ni  is  Poisson  distributed  with  mean 
and  variance  of  A^  +  Ab 


The  Raman  spectroscopy  sensor  system  transmits  a  laser 
pulse  and  produces  a  measured  Raman  spectrum  from  the 
energy  scattered  by  the  chemical  substance.  The  spectrum 
is  spread  across  the  bins  of  a  CCD  detector.  The  response 
on  each  bin  corresponds  to  the  amount  of  energy  scattered 
at  a  particular  frequency  or  wave  number. 

Let  y  e  Tv M  denote  a  measured  spectrum  with  values  at 
Mbins 


S  ,  b 

n,  +n; 


^Poisson  (ni  +»,b;^  +^b)- 


The  vector  measurement  model  is 


y=n  +n  +g, 


(9) 

(10) 


where  (ns,nb,g)e  are  defined  similarly  as  in  (1). 


Large  Signal  Approximation 


y  :=  [p,  v2  •••  yMl  (!) 

where  is  used  to  define  a  quantity.  The  measurement 
model  [16-17],  [11]  for  each  element  ofy  is  described  by 

yt  =  n-  +  n,b  +  gt ,  i  =  1,2, (2) 


If^  is  large,  then  iy  +  ly  is  well  approximated  by 
Gaussian  distributions 


n*  +  n'>~N{n*  +  n';-,A?i+Abi,Xi+A)).  (11) 
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®  =  PA. 


Using  (2)  and  (11) 

y,  ~  N( n°  +  n) b  +  gi;A)  +  X)  +  +  X)  +  a2).  (12) 

Alternatively, 

yi=^+^  +  m  +  v„  (13) 

vt~N(  0,^+^  +  a2).  (14) 

Using  the  large  signal  approximation,  (13)  can  be  written 
as 

y  =  +  >^b  +  m  +  v,  (15) 

where 

m:=m[l  1,...,1]',  (16) 

V:=[%  (17) 

Ab  :=[/lb  (18) 

v  ~  N(\;  0Mx1,R),  (19) 


(28) 

Not  all  photons  that  hit  the  CCD  array  are  converted  to 
photoelectrons.  The  quantum  efficiency  or  flat-field 
response  varies  along  the  CCD  array.  This  non-uniform 
detector  efficiency  is  modeled  by 

^:=A(H  (29) 

where  /?;is  known  from  calibration  measurements.  We 
can  write  (29)  in  the  matrix  form 

r=cx,  (30) 

where 

C:=BC*  =  BPA,  (31) 

B  =  diag(/?,, (32) 

Under  the  large  signal  approximation,  substitution  of  (30) 
in  (15)  gives 

y  =  Cx  +  lb+m  +  v.  (33) 


R diag(/l|  +  +  a~ , . . . ,  XM  +  +  (T~ ).  (20) 

3  Measurement  Function  for 
Raman  Spectra 

Suppose  we  have  A  reference  spectra  {s ;  e  M  in  our 

library  corresponding  to  A  chemical  substances.  Then  the 
true  target  spectra  s  can  be  expressed  as  a  linear 
combination  of  the  reference  spectra  by 

S  =  Z  xj*r  (21) 

7=1 

We  can  write  (21)  in  the  matrix  form 

s  =  Ax,  (22) 

where 

x .— [Xi  x2  ...  xN] ,  (23) 

Xj>  0,  j  =1,2,..., N,  (24) 

A  :=  [Sj  s2  ...  sN].  (25) 

Then 

V  =  PS,  (26) 


Define  the  new  measurement  vector  z 

^  b 

z  :=  y  —  k  -  m. 


Then 


z  =  C  x  +  v. 


(34) 

(35) 


Thus,  under  the  large  signal  approximation,  the 
measurement  model  is  linear  with  additive  Gaussian 
measurement  noise.  An  estimate  x  of  x  can  be  obtained 
using  the  maximum  likelihood  estimator  (MLE)  [4],  [10] 
or  weighted  least  squares  (WLS)  [4],  [10]  with  the 
nonnegativity  constraint  (24).  Thus,  the  estimation 
problem  is  a  constrained  estimation  problem  due  to  (24). 
Use  of  classical  MLE  or  WLS  would  yield  approximate 
results. 


4  Raman  Spectra  Estimation  Algorithms 

This  paper  addresses  estimation  of  Raman  spectra  under 
the  large  signal  assumption.  Future  work  will  address  the 
more  general  case  where  the  large  signal  assumption  is 
not  valid. 


Since  the  measurement  model  in  (35)  is  linear  with 
additive  Gaussian  noise,  the  estimates  from  the  MLE  and 
WLS  are  the  same  provided  the  weight  matrix  W  in  WLS 
is  equal  to  R_1  [4],  [10].  Then  using  (20), 


where  P  is  the  M  X  M  point  spread  function  matrix  of  the  W  :=  diag(  vv, ,  w2 , . . . ,  wM  ),  (36) 

diffraction  grating  used  to  spread  the  spectral  energy 

across  the  CCD’s  bins.  Substitution  of  (22)  in  (26)  gives  wi  :=  1  /(A)  +  /lb  +  CT2),  i  =  1,2,..., M.  (37) 


V  =  <Px, 


(27)  The  true  weights  in  (37)  are  not  known  and  must  be 
estimated.  Our  future  work  will  address  this  issue. 


where 
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The  cost  function  for  the  parameter  estimation  problem  is 

J(x)  :=  (z  -  C  x)'W(z  -  C  x).  (38) 

We  can  rewrite  (38)  as 

J(x):=(q-Dx)'(q-Dx),  (39) 

where  the  weighted  measurement  vector  r|  e  9t"  and 
weighted  measurement  matrix  De  9v  MxN  are  defined  by 

r/i  :=  i  =  1,2,..., M,  (40) 

dy^4^icu>  *  =  1,2,..., Af,  y  =1,2,..., AT.  (41) 

4.1  Classical  Weighted  Least  Squares 

Classical  WLS  (CWLS)  [4],  [10]  solves  the  following 
problem  without  the  nonnegative  constraint  (24) 

minJ(x).  (42) 

x  v  7 

Since  the  CWLS  does  not  enforce  the  nonnegative 
constraint,  the  estimate  obtained  using  the  CWLS  is 
expected  to  have  lower  accuracy  compared  to  the 
nonnegative  WLS  (NNWLS). 

4.2  Nonnegative  Weighted  Least  Squares 

For  the  current  measurement  model  (35),  the  NNWLS  or 
nonnegative  MLE  (NNMLE)  solves  the  problem 

minJ(x).  (43) 

x>0 

There  are  three  commonly  used  algorithms  for  solving  the 
NNWLS  problem  in  (43).  These  algorithms  were 
proposed  by  Lawson  and  Hanson  [9]  (LHNNLS),  Bro  and 
De  Jong  [1],  and  Van  Benthem  and  Keenan  [18] 
(VKNNLS).  Although  these  three  algorithms  considered 
the  same  measurement  variances  for  all  measurements,  i.e. 
W  =  cr2I ,  they  can  be  easily  modified  to  handle  non- 

uniform  weights.  The  latter  two  algorithms  are 
improvements  over  the  original  algorithm  of  Lawson  and 
Hanson  for  handling  multiple  measurement  vectors, 
zk  e  9LM,  k  =  1,2, It  can  also  be  shown  that  for 

multiple  measurement  vectors  {zk  G  )RM  }  ,  solving  the 

problem  by  processing  one  measurement  vector  at  a  time 
is  much  less  efficient  than  collecting  a  number  of 
measurement  vectors  and  processing  them  at  once  using  a 
column  parallel  algorithm  [18].  Next,  we  summarize  the 
three  NNWLS  algorithms. 


4.2.1  Lawson-Hanson  Algorithm 

The  standard  algorithm  for  computing  XNNWLS  is  that  of 

Lawson  and  Hanson  [9],  which  is  included  in  the 
VIATLAB  function,  "lsqnonneg.” 

Algorithm  1  LHNNWLS  [9] 

Let  Sp  and  Sa  denote  the  passive  and  active  index  sets, 
respectively. 

Given  D,  q,  find  x  that  solves  (43)  with  (39). 

1.  Set  S^NULL,  Sa  =  {1,2  ...V},  and  x  =  0. 

2.  Compute  the  V-vector  w  =  D’  (q  -  Dx). 

3.  If  the  set  S&  is  empty  or  if  w  .  <  0  for  all  j  e  Sa , 

then  go  to  Step  12. 

4.  Find  an  index  t  e  5a  such  that 

w,  =  maxjw j  :  j  e  Sa} 

5.  Move  the  index  t  from  set  S  to  set  S  ■ 

a  p 

6.  Let  Dp  denote  the  m2xN  matrix  defined  by 


Column  j  of  D  _  Jcoiumn  j  of  D  if  y  e  5p  | 

P  1°  J' 

Compute  the  N- vector  q ,  as  a  solution  of  the  CLS 
problem  Dpx  s  q. 

Note  that  only  the  components  x j,j  e  Sp  ,  are 
determined  by  this  problem. 

Define  x,  =  0  for  ye  A. 

J  J  a 

7.  If  x j  >  0  for  all  j  e  5p ,  set  x  =  q  and  go  to  Step  2. 

8.  Find  an  index 

qe  Sp  such  that  xq /(xq  -r)q)  = 
min{  Xj  / (Xj  -Xj):  Xj  <  0,  j  e  Sp  }. 

9- Set  a  =  xq/(xq-xq). 

10.  Set  x  =  x  +  «(q-x). 

1 1 .  Move  from  set  Sp  to  set  Sa  all  indices  j  e  Sp  for 


which  Xj  =  0  .  Go  to  step  6. 

12.  Computation  completed 


XNNWls=  x  *s  solution  such  that: 


j  Xj  >  0,/e  Sp] 
[xj=0,jeSt  J' 


End  Algorithm  1  LHNNLS 


This  algorithm  has  some  disadvantages  for  multiple 
measurement  vectors  because  of  the  repeated  solution  of 
the  CLS  problem  Dpx  s  q  in  Step  6  within  each  of  the 

iterations. 


2242 


4.2.2  Fast  Combinatorial  NNWLS 
Algorithms  (FCNNWLS) 


The  efficiency  of  the  NNWLS  can  be  improved  by 
processing  multiple  measurement  vectors  in  a  block, 
(zj,z2,...,zt) ,  for  unweighted  measurements,  and 


(Ill. 111. 

x\k),  for  weighted  measurements. 

For  multiple 

measurement  vectors,  define  the  weighted  measurement 

matrix  G  and  parameter  matrix  X 

G  :=  [iii  *12  ■••%], 

(44) 

X  :=  [x,  x2...xj. 

(45) 

Then  the 

NNWLS  problem  for  multiple 

measurement 

vectors  is 

min  DX-G  2 , 

x>o 

(46) 

where  F  in  (46)  represents  the  Frobenius  norm  [9].  This  is 
possible  since  common  least  squares  (LS)  problems  across 
multiple  observation  vectors  can  be  processed 
simultaneously.  The  resulting  modification  to  LF1NNLS  is 
referred  to  as  the  column  parallel  form  of  processing  the 
passive  set  associated  with  multiple  measurement  vectors. 
The  redundant  computations  can  be  taken  advantage  of 
and  the  pseudo  inverse  [9]  of  D;,  in  Step  6  can  be 

computed  for  some  of  the  multiple  measurement  vectors 
with  common  structure  [6,7,18].  All  measurement  vectors 
that  correspond  to  a  common  pseudoinverse  are  grouped 
together.  Thus,  the  pseudoinverse  is  computed  only  once 
for  each  group  of  measurement  vectors  sharing  this 
common  pseudoinverse. 

Modifications  to  algorithm  LF1NNLS  are  made  to  the  CLS 
part  of  the  code  in  Step  6.  These  entail  finding  the 
multiple  measurement  vectors  with  a  common 
pseudoinverse  and  solving  the  CLS  problem.  This 
pseudoinverse  may  be  used  in  subsequent  computations  as 
well.  As  an  example,  given  3  measurement  vectors,  the 
number  of  pseudoinverse  computations  may  be  reduced 
from  7  to  4  (see  [20]  for  details).  Two  pseudoinverse 
computations  are  common  to  two  of  the  three  iterations 
and  one  of  those  is  used  twice  in  the  second  iteration.  The 
sequence  is 

1.  {1,1,1} 

2.  {2,3,2} 

3.  {2,3,4} 

Where  the  same  pseudoinverse  is  applied  to  all  three 
columns  in  step  1,  one  less  pseudoinverse  is  required  in 
step  2,  and  two  of  those  are  reused  in  step  3.  The  savings 
in  computations  will  become  even  more  significant  as  the 
number  of  measurement  vectors  becomes  larger. 


4.2.3  Block  Pivoting  NNWLS  (BPNNWLS) 

Further  computational  savings  can  be  realized  by 
generalizing  the  active  set  method  of  LF1NNLS,  which  is  a 
single  principal  pivoting  algorithm,  to  a  block  principal 
pivoting  algorithm  [8].  This  method  modifies  LF1NNLS 
using  a  block  exchange  rule  to  move  blocks  of  columns 
from  the  “passive”  to  “active”  sets.  The  quotes  are  due  to 
the  fact  that  the  sets  employed  in  the  block  principal 
pivoting  method  do  not  necessarily  correspond  to  the 
active  and  passive  sets  of  algorithm  LF1NNLS. 

4.2.4  Weighted  Generalized  Likelihood 
Ratio  Test  (WGLRT)  Algorithm 

The  subspace  version  of  the  generalized  likelihood  ratio 
test  (GLRT)  [5]  has  been  applied  to  analyze  Raman 
spectroscopy  data  [14-15].  In  this  formulation,  the 
likelihood  function  p( z;  x,  R,  H] )  for  the  measurement 
model  (35)  is  calculated  under  the  hypothesis  Hx ,  where 
all  reference  spectra  are  included  in  constructing  C  and  x 
is  the  NNWLS  estimate.  Then  the  zth  reference  spectra  is 
removed  and  the  likelihood  function  p(z;  x0,  R0,  H0)for 
the  measurement  model 


H  '  z  —  C  x  +  v 

11  0  •  ^0  A0  ~  Y  0’ 

(49) 

v0~A(v0;0,R0), 

(50) 

is  calculated  under  the  hypothesis  H0  ,  where  C0  is  the 

corresponding  measurement  matrix  and  x0  is  the 

NNWLS  estimate.  If  the  log-likelihood  ratio  exceeds  a 
threshold,  i.e., 

In  p( z;  X,  R,  Hx )  -  In p(z;  x0 ,  R0 ,  H0 )  >  a,  (51) 

then  the  zth  substance  is  assumed  to  contribute  to  the 
target  substance.  This  process  is  repeated  for  each 
reference  spectra  and  the  indices  {il,i2,...,ip}  are 
determined  for  which  the  test  (51)  succeeds.  Using  these 
indices,  the  measurement  matrix  C  is  formed  and  the 

NNWLS  estimate  x  is  calculated. 

p 

5  Numerical  Simulation  and  Results 

We  used  67  reference  Raman  spectra,  {s;.  e  9tM}^=1 , 

N=61.  All  the  components  of  x  were  zero  except 
jc60  =1.0.  Each  spectrum  has  values  at  M=  1024  bins.  In 

the  Monte  Carlo  simulations,  the  mean  and  variance  of  the 
Gaussian  measurement  noise  are  10  and  225,  respectively. 
We  used  a  constant  value  of  256  for  the  Poisson  parameter 
A’’  for  all  bin  values.  We  then  calculated  Asby  selecting 

and  substituting  a  true  x  vector  into  (30).  Figure  1  shows 
the  variation  of  the  measurement  error  variance  with  bin 
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i  =  1,2, 


(53) 


index  for  the  current  scenario.  We  observe  that  the 
measurement  error  variance  changes  significantly  with  the 
bin  index. 


Figure  1.  Variation  of  measurement  error  variance  with 
bin  index. 

One  hundred  Monte  Carlo  trials  were  used  to  calculate 
measures  of  performance  (MoP)  for  each  spectral 
estimation  algorithm.  The  equations  and  resulting  MoP 
results  are  presented  here. 


Zm,i  '  Ztn,i  Zm,i'> 


The  bias  error  of  the  measurement  residual  at  the  zth  bin 
and  the  overall  bias  of  the  measurement  residual  are 
defined,  respectively,  by 


Kr-=irl 

M,  ^ 


J  m,i 9 


(54) 


b 


z 


1 

MMS 


M  M, 

X2X,o 


i=l  m= 1 


(55) 


The  RMSE  of  the  measurement  residual  at  the  /"'bin  and 
the  overall  RMSE  for  the  measurement  residual  are 
defined,  respectively,  by 


RMSE  ,  := 


1 

M, 


£ 


i  =  (56) 


Let  Ms  be  the  total  number  of  Monte  Carlo  simulations 
and  x  the  estimate  of  x  in  the  /71th  Monte  Carlo 

m 

simulation.  The  estimation  error  in  the  yth  component  of 
x  in  the  771th  Monte  Carlo  simulation  is  defined  by 

x  :=  x  —  x  ,  i  =  l,2,...,N.  (47) 

m,j  7  m,j  ’  J  j-5^-5 . . .  9  j.  »  .  \  / 


RMSEr 


1 

MMS 


M  Ms 

I  2Xi 


i=l  m= 1 


1/2 


(57) 


Table  1  summarizes  unweighted  bias  and  RMS  error 
results  for  both  parameter  estimation  and  measurement 
residual  using  100  Monte  Carlo  trials  and  one 
measurement  vector.  Recall  that  unweighted  refers  to  the 
use  of  a  weight  matrix  equal  to  the  identity  matrix. 


The  bias  error  for  the  jth  coefficient  and  overall  bias  error 
for  the  coefficients  are  defined,  respectively,  by 


bx  ■  :=  — 

x’]  M. 


2> 


m,j ? 


j  =  \,2,...,N, 


(48) 


b 


X 


1 

mTs 


N  Ms 

ZIXi- 


7=1  m= 1 


(49) 


The  root  mean  square  error  (RMSE)  for  the  jth  coefficient 
and  the  overall  RMSE  for  the  coefficients  are  defined, 
respectively,  by 


RMSE,,.  := 


m— 1 


1/2 


j  =  1,2, (50) 


RMSE,  := 


N 


NM, 


s  7=1  m= 1 


(51) 


Let  Zm  .  denote  the  predicted  measurement  at  the  /"'  bin  in 
the  777  th  Monte  Carlo  simulation.  Then 


4,7  :=(CiJ„  i  =  1,2, (52) 

The  measurement  residual  at  the  /"’bin  in  the  /77th Monte 
Carlo  simulation  is  defined  by 


Table  1.  Overall  errors  for  unweighted  versions 
of  the  algorithms. 


Algorithm 

Estimation 

Residual 

Bias 

RMSE 

Bias 

RMSE 

CLS 

-1.5904e-5 

0.0928 

-0.1041 

7.4593 

NNLS 

1.4883e-4 

0.0031 

0.9737 

3.1252 

FCNNLS 

1.4883e-4 

0.0031 

0.9737 

3.1252 

BPNNLS 

1.4883e-4 

0.0031 

0.9737 

3.1252 

NNGLRT 

8.8715e-13 

4.6835e-4 

-1.04e-6 

0.7288 

Table  2  summarizes  the  bias  and  RMS  error  results  when 
the  algorithms  use  the  weight  matrix  defined  in  (36-37). 
Comparing  the  results  between  Table  1  and  Table  2  shows 
a  significant  reduction  in  the  measurement  residual  bias 
error  and  RMSE  when  the  weight  matrix  is  used. 


Table  2.  Overall  errors  for  weighted  versions 
_ of  the  algorithms. _ 


Algorithm 

Estimation 

Residual 

Bias 

RMSE 

Bias 

RMSE 

CWLS 

-1.6241e-5 

0.0918 

-0.0053 

0.8240 

NNWLS 

1.1793e-4 

0.0028 

0.0421 

0.1801 

FCNNWLS 

1.1793e-4 

0.0028 

0.0421 

0.1801 

BPNNWLS 

1.1793e-4 

0.0028 

0.0421 

0.1801 

NNWGLRT 

1.102e-14 

4.322e-4 

0.0018 

0.0256 

The  results  for  parameter  estimation  are  nearly  the  same 
for  the  first  four  algorithms.  In  addition,  the  NNWGLRT 
algorithm  has  the  best  performance. 
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Tables  3  and  4  show  the  average  CPU  times  for  100 
Monte  Carlo  trials.  The  LS-based  methods  have  lower 
CPU  times  than  the  likelihood-based  method.  For 
multiple  measurement  vectors  the  block  pivoting  method 
outperformed  the  other  constrained  least  squares  methods. 
The  CLS  method  is  faster  but  the  estimation  and  residual 
errors  are  very  large  compared  to  the  constrained  methods. 

Table  3.  Average  CPU  time  over  100  Monte  Carlo  trials. 


Algorithm 

CPU  Time  (seconds) 

Unweighted 

Weighted 

CLS 

3.18 

3.26 

NNLS 

3.20 

3.23 

FCNNLS 

3.76 

3.80 

BPNNLS 

3.32 

3.34 

NNGLRT 

11.68 

11.77 

Table  4.  Average  CPU  time  over  100  Monte  Carlo  trials. 
Algorithms  used  to  process  multiple  measurement  vectors. 


Algorithm 

CPU  Time  (seconds) 

Unweighted 

Weighted 

20  MVs 

20  MVs 

100  MVs 

CLS 

2.64 

2.69 

2.84 

NNLS 

4.89 

5.07 

11.08 

FCNNLS 

5.54 

5.67 

13.02 

BPNNLS 

3.31 

3.38 

5.90 

These  results  demonstrate  that  great  care  needs  to  be  taken 
when  specifying  the  measurement  model  since  this 
impacts  the  formulation  of  the  detection  algorithms.  As 
can  be  seen  from  the  tables  and  the  graphs,  the  residual 
errors  can  be  reduced  significantly  with  the  proper 
measurement  model  and  measurement  error  covariance. 
Note  the  decrease  in  the  residual  error  bias  and  RMSE  in 
Tables  1  and  2.  Plotting  MoP  values  versus  bin  index 
averaged  over  the  100  Monte  Carlo  trials  in  Figures  2,  3, 
and  4  clearly  illustrates  the  decrease. 

Figure  2  compares  the  bias  errors  of  measurement 
residuals  for  the  CLWS,  NNWLS,  and  NNWGLRT 
algorithms.  The  NNWGLRT  algorithm  has  nearly  zero 
bias. 


bin  number 

Figure  2.  CWLS,  NNWLS  and  NNWGLRT  residual  error  bias. 


performance  with  respect  to  residual  error  RMSE.  The 
NNWGLRT  results  exhibit  almost  no  residual  error.  It  is 
also  interesting  to  note  that  the  unweighted  NNGLRT 
algorithm  performs  better  than  the  NNLS. 


Figure  3.  CLS,  CWLS,  NNLS,  and  NNWLS  residual  error  RMSE. 


bin  number  (100  to  350) 

Figure  4.  NNLS,  NNWLS,  NNGLRT,  and  NNWGLRT  residual  error 
RMSE  zoomed  to  bins  100  to  350. 


Figures  5  and  6  show  the  bias  and  RMSEs  for  parameter 
estimation  using  the  CWLS,  NNWLS,  and  NNWGLRT 
algorithms.  We  observe  that  the  use  of  nonnegative 
constraints  significantly  improves  the  accuracies  of 
parameter  estimation. 
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Figure  5.  Estimation  error  bias. 


mixing  coefficient 


Figure  6.  Parameter  estimation  RMSE. 


In  Figure  3,  we  observe  that  CLS  has  the  highest  RMSE 
for  measurement  residual  and  NNWLS  has  the  lowest. 

We  also  note  that  CWLS  outperforms  NNLS.  This 
observation  underscores  the  importance  of  including  error 
sources  in  the  measurement  model.  Figure  4  shows  that 
both  the  NNWLS  and  NNWGLRT  exhibit  good 


6  Conclusions 

This  paper  compared  five  supervised  algorithms  for 
estimating  Raman  spectra  in  the  large  signal  domain.  In 
this  domain,  the  measurement  model  can  be  approximated 
as  being  linear  with  additive  Gaussian  noise.  The 
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measurement  error  variances  vary  significantly  with  bin 
index.  The  parameter  estimation  problem  is  constrained 
because  the  elements  of  the  parameter  vector  must  be 
nonnegative.  We  have  presented  numerical  results  for  the 
bias  error  and  RMSE  of  the  estimated  parameter  and 
measurement  residual. 

Since  the  measurement  error  variances  vary  significantly 
with  bin  index,  it  is  important  to  use  the  correct  weights  or 
measurement  error  variances  in  parameter  estimation.  It  is 
also  important  to  enforce  the  nonnegativity  constraint  in 
the  estimation  of  the  mixing  coefficients.  The  NNWLS, 
FCNNWLS,  and  BPNNWLS  algorithms  yield  nearly  the 
same  accuracy  and  NNWGLRT  has  the  best  accuracy,  but 
the  worst  computational  performance.  The  CPU  times  of 
NNWLS,  FCNNWLS,  and  BPNNWLS  are  also  similar 
when  only  one  measurement  vector  is  processed  at  a  time. 
When  a  block  of  data  is  processed  together,  BPNNWLS 
shows  considerable  computational  advantage. 

Our  future  work  will  focus  on  three  areas.  First,  the  more 
general  case  when  the  large  signal  approximation  is  not 
valid  will  be  investigated.  Second,  a  more 
computationally  efficient  NNGLRT  algorithm  will  be 
developed.  Third,  a  similar  investigation  and  comparison 
with  unsupervised  algorithms  will  be  carried  out  so  that 
their  potential  for  Chemical/Biological  detection 
capability  can  be  evaluated. 
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